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Abstract 

An  approach  based  on  a  discontinuous  Galerkin  discretization  is  proposed 
for  the  Bhatnagar-Gross-Krook  model  kinetic  equation.  This  approach  al¬ 
lows  for  a  high  order  polynomial  approximation  of  molecular  velocity  dis¬ 
tribution  function  both  in  spatial  and  velocity  variables.  It  is  applied  to 
model  normal  shock  wave  and  heat  transfer  problems.  Gonvergence  of  solu¬ 
tions  with  respect  to  the  number  of  spatial  cells  and  velocity  bins  is  studied, 
with  the  degree  of  polynomial  approximation  ranging  from  zero  to  four  in 
the  physical  space  variable  and  from  zero  to  eight  in  the  velocity  variable. 
This  approach  is  found  to  conserve  mass,  momentum  and  energy  when  high 
degree  polynomial  approximations  are  used  in  the  velocity  space.  For  the 
shock  wave  problem,  the  solution  is  shown  to  exhibit  accelerated  conver¬ 
gence  with  respect  to  the  velocity  variable.  Gonvergence  with  respect  to  the 
spatial  variable  is  in  agreement  with  the  order  of  the  polynomial  approxima¬ 
tion  used.  For  the  heat  transfer  problem,  it  was  observed  that  convergence  of 
solutions  obtained  by  high  degree  polynomial  approximations  is  only  second 
order  with  respect  to  the  resolution  in  the  spatial  variable.  This  is  attributed 
to  the  temperature  jump  at  the  wall  in  the  solutions.  The  shock  wave  and 
heat  transfer  solutions  are  in  excellent  agreement  with  the  solutions  obtained 
by  a  conservative  finite  volume  scheme. 


*  Corresponding  author.  Phone:  -1-1  818-677-2645.  Fax:  -|-1  818-677-3634 
Email  addresses:  alexander.alekseenko@csun.edu  (Alexander  Alekseenko  ), 
ngimel@gmail.com  (Natalia  Gimelshein) ,  gimelshe@usc.edu  (Sergey  Gimelshein) 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


1  INTRODUCTION 


Key  words: 

high  order  methods,  discontinuous  Galerkin  methods,  model  kinetic 

equations,  shock  wave,  heat  transfer 

PACS:  47.45.Ab,  02.60.-x,  47.61. Cb,  44.05.+e 


1.  Introduction 

Over  the  last  decade,  numerical  modeling  of  flows  in  gas-driven  micro- 
and  nanoscale  devices  has  drawn  much  attention  due  to  its  application  to 
sensors,  actuators,  hlters,  pumps,  flow  control  systems,  and  so  forth.  For 
many  of  these  devices,  development  stems  from  trial  and  error  approaches 
htted  around  device  fabrication.  For  their  development  to  take  a  leap  for¬ 
ward,  accurate  and  efficient  numerical  modeling  is  necessary  in  areas  where 
flow  diagnostics  may  be  limited  or  impossible.  The  choice  and  use  of  the 
numerical  approach  suitable  to  microscale  flows  is  primarily  related  to  a 
number  of  fluid  dynamic  effects  peculiar  to  these  flows.  These  effects,  such 
as  velocity  slip,  temperature  jump,  thermal  creep,  and  viscous  heating,  dom¬ 
inate  in  many  microdevices  based  on  microelectromechanical  technologies. 
Due  to  the  small  characteristic  dimensions  of  microscale  flows — typically  on 
the  order  of  tenths  of  micrometers  to  millimeters — heat  conductivity  and 
diffusivity  play  a  very  important  role.  Large  surface-to-volume  ratios  imply 
that  gas-surface  interactions  are  very  important.  The  gas  mean  free  paths 
are  comparable  to  characteristic  flow  dimensions,  which  results  in  signihcant 
deviation  from  equilibrium.  This  means  that  in  many  cases  conventional 
computational  fluid  dynamics  methods  such  as  the  solution  of  Navier-Stokes 
equations  are  not  applicable.  Instead  modeling  must  be  based  on  kinetic  gas 
theory. 

Kinetic  gas  theory  describes  gas  properties  through  the  distribution  func¬ 
tion  of  molecular  velocities.  The  governing  equation  for  the  velocity  distri¬ 
bution  function  is  the  Boltzmann  equation.  It  expresses  the  variation  of  the 
distribution  function  due  to  molecular  free  flight,  action  of  external  forces  and 
intermolecular  collisions.  For  a  single  species  monatomic  gas,  the  Boltzmann 
equation  has  the  form 

f  +  "'i  +  -  f  h'iSbdbdednu  (1) 

where  the  velocity  distribution  function  /  is  dehned  by  the  condition  that 
f(t,x,u)  dxdu  is  the  number  of  molecules  at  time  t  with  velocities  between 
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u  and  u+  du  and  coordinates  between  x  and  x+  dx.  Here,  X  is  the  external 
force,  g  =  \u  —  ui\  is  the  magnitnde  of  relative  velocity  of  colliding  molecules, 
and  b  and  e  are  geometric  impact  parameters. 

The  Boltzmann  equation  is  a  nonlinear  integro-differential  equation  amenable 
to  analytical  solution  only  for  a  small  number  of  special  cases  of  collisionless 
or  spatially  homogeneous  problems.  Difficulties  encountered  in  numerical 
solution  of  the  Boltzmann  equation  are  primarily  attributed  to  the  multidi¬ 
mensional  phase  space — physical  coordinates  and  velocity  coordinates — and 
the  multidimensional  collision  integral  in  the  right-hand  side.  For  three- 
dimensional  flows  of  monatomic  gases  the  collision  integral  involves  integra¬ 
tion  over  a  hve-dimensional  domain.  This  limits  the  direct  numerical  inte¬ 
gration  of  the  Boltzmann  equation  [1]  to  simplihed  collision  models,  such 
as  hard  sphere  molecules,  and  one-  and  two-dimensional  flow  problems.  The 
molecular  dynamics  method  [2]  is  also  applicable  to  the  solution  of  the  Boltz¬ 
mann  equation,  but  it  is  even  more  limited  in  terms  of  flow  dimensionality 
and  problem  size.  The  most  powerful  and  widely  used  approach  to  the  so¬ 
lution  of  the  Boltzmann  equation  is  currently  the  direct  simulation  Monte 
Carlo  (DSMC)  method  [3]. 

Although  the  DSMC  method  is  more  efficient  than  the  two  former  ap¬ 
proaches,  it  still  suffers  from  high  computational  cost  compared  to  conven¬ 
tional  continuum  CFD  methods.  The  computational  cost  limitation  of  the 
DSMC  method  are  especially  severe  at  low  Knudsen  numbers;  in  many  cases 
it  makes  this  method  an  unacceptable  choice.  This  is  the  case  for  three- 
dimensional  low  speed  flows  where  the  computational  cost  is  impacted  by 
flow  dimensionality,  the  long  time  to  reach  steady  state,  the  low  signal-to- 
noise  ratio,  multiple  physical  scales  and  other  factors. 

There  are  several  alternative  DSMC-based  approaches  proposed  to  deal 
with  the  problem  of  low  signal-to-noise  ratio  that  allow  reduction  in  macropa¬ 
rameter  sampling  time  compared  to  the  standard  DSMC  method  [4,  5,  6,  7]. 
Although  all  these  techniques  signihcantly  reduce  steady-state  time  averaging 
cost,  they  do  not  deal  with  the  computational  cost  associated  with  long  times 
to  reach  steady  state  typical  for  low  speed  micro-  and  nanoscale  flows.  This 
time  is  usually  comparable  with  or  larger  than  the  time  needed  to  sample 
macroparameters,  and  thus  any  attempt  to  shorten  the  latter  is  bounded  by 
the  former.  In  addition,  these  methods  do  not  address  the  problem  of  particle 
correlations  and  related  solution  accuracy,  peculiar  for  particle  methods  and 
greatly  amplihed  for  low-speed  flows.  Finally,  particle  approaches  are  poorly 
suited  for  unsteady  fluid-thermal  coupling  due  to  extremely  high  statistical 
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scatter  of  instantaneous  ensemble-averaged  heat  fluxes  from  gas  to  surface. 

This  means  that  current  applications  of  particle-based,  statistical  ap¬ 
proaches  to  the  solution  of  the  Boltzmann  equation  for  modeling  low-speed 
microscale  gas  flows  are  fairly  limited.  A  plausible  numerical  alternative  is  a 
deterministic  solution  of  several  available  simplihed  forms  of  the  Boltzmann 
equation,  known  as  model  kinetic  equations.  Two  of  the  most  known  model 
kinetic  equations,  the  Bhatnagar-Gross-Krook  (BGK)  [8]  and  the  ellipsoidal 
statistical  (ES)  [9]  kinetic  models,  use  a  non-linear  relaxation  term  instead 
of  the  full  Boltzmann  collision  integral.  In  spite  of  the  simpler  collision  term, 
both  models  possess  the  same  collision  invariants  as  the  Boltzmann  equation. 

The  main  objective  of  this  work  is  to  present  a  new  approach  to  the  so¬ 
lution  of  the  BGK  model  kinetic  equation  based  on  a  Galerkin  discretization 
of  both  physical  and  velocity  space  and  to  analyze  its  applicability  to  model 
microscale  gas  flows.  Advantages  and  disadvantages  of  the  discontinuous 
Galerkin  (DG)  method  of  solution  of  model  kinetic  equations  are  discussed. 
For  two  classical  problems  of  rarehed  gas  dynamics — normal  shock  wave  and 
heat  transfer  between  parallel  plates — the  DG  solutions  are  compared  with 
the  solutions  obtained  by  a  hnite  volume  (FV)  approach. 

2.  Challenges  of  the  numerical  solution  of  the  BGK  equation 

The  BGK  equation  is  a  kinetic  model  obtained  from  an  intuitive  simpli- 
hcation  of  the  collision  integral  in  the  Boltzmann  equation.  Similar  to  the 
Boltzmann  equation,  it  is  solved  for  the  multi-dimensional  molecular  velocity 
distribution  function.  In  the  BGK  model  the  collision  operator  is  replaced 
with  a  simpler  operator 

«^(/o  -/), 

where  z/(t,  x)  is  the  collision  frequency  and  fo(t,  x)  is  the  Maxwellian  equilib¬ 
rium  distribution  function 

/o  =  n(t,f)(27ri?T(t,T))~^/^exp(^-^^  2RT^^^  )' 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


2  CHALLENGES  OF  THE  NUMERICAL  SOLUTION  OF  THE  BGK 

EQUATION 


The  gas  density,  n,  bulk  velocity,  u  and  temperature,  T,  are  defined  as  follows: 

n{t,x)=  /  f{t,x,u)du; 

2r3 

n{t,x)u{t,x)  =  /  uf{t,x,u)du] 

JRS 

n{t,x)T{t,x)  =  ^  f  {u-uff{t,x,u)du,  (3) 

6K  J^3 

where  R  is  the  gas  constant.  The  BGK  model  satisfies  the  mass,  momentum, 
and  energy  conservation  laws  and  the  Boltzmann’s  H-theorem  expressing  the 
increase  of  entropy  of  the  gas  under  consideration.  Kinetic  models  should 
also  reproduce  the  gas  transport  coefficients — viscosity,  thermal  conductivity 
and  species  diffusivity — resulting  from  the  Boltzmann  equation.  The  primary 
advantage  of  this  numerical  alternative  is  its  relatively  high  computational 
efficiency.  Previous  application  of  the  solution  of  the  model  kinetic  equations 
to  low-speed  microscale  gas  flows  showed  an  improvement  in  numerical  effi¬ 
ciency  of  more  than  two  orders  of  magnitude  using  a  deterministic  solution 
of  model  kinetic  equations  was  used  instead  of  the  DSMC  method  [10,  11]. 

Current  approaches  to  the  solution  of  model  kinetic  equations  include 
either  finite  difference  [10]  or,  much  more  frequently,  finite  volume  methods 
[12,  13,  14].  Finite  volume  approaches  are  generally  more  robust  and  may  be 
applied  to  complex  two-  and  three-dimensional  geometries.  There  is  also  a 
methodology  that  enforces  conservation  laws  [12].  However,  there  are  several 
problems  related  to  the  application  of  finite  volume  approaches  to  model 
low-speed  microfiows.  First,  their  efficiency,  while  much  higher  than  that  of 
DSMC,  is  still  much  lower  than  the  efficiency  of  continuum  CFD  solvers.  This 
is  primarily  related  to  the  multi-dimensionality  of  the  velocity  distribution 
function  which  depends  on  three  spatial  coordinates  and  three  components 
of  velocity.  The  existing  finite  volume  solvers  for  model  kinetic  equations 
are  typically  second  order  in  physical  space.  Their  approximation  in  velocity 
space  is  generally  based  on  piece-wise  constant  functions,  also  known  as  the 
discrete  ordinate  approach.  As  a  result,  a  large  number  of  velocity  bins 
and  spatial  cells  is  usually  needed  to  obtain  a  converged  solution.  This 
significantly  limits  the  application  of  these  solvers  to  three-dimensional  flows 
with  a  considerable  degree  of  non-equilibrium. 

The  second  problem  is  related  to  virtual  cells  behind  the  surface  that  are 
necessary  in  the  implementation  of  finite  volume  schemes.  The  use  of  such 
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cells  impacts  the  accuracy  of  calculation  of  fluxes  at  gas-solid  interfaces. 
Errors  caused  by  the  use  of  virtual  cells  not  only  affect  gas  properties  near 
the  surface,  but  also  propagate  into  the  gas  volume.  The  property  affected 
most  is  the  gas  velocity.  An  example  of  undesirable  numerical  effects  induced 
by  virtual  cells  is  given  in  Figure  1,  where  the  x  velocity  and  streamlines  are 
shown  around  a  radiometer  vane  immersed  in  a  6  Pa  argon  gas  and  heated 
to  415  K  on  the  left  side  and  395  K  on  the  right  side.  The  radiometer  is 
centered  at  (0, 0).  The  lower  boundary  is  the  symmetry  plane,  and  the  other 
boundaries  are  chamber  walls  kept  at  300  K.  It  is  seen  that  near  the  chamber 
walls  there  are  hnite  (about  0.1  m/s)  flow  velocities  in  the  direction  normal 
to  the  wall.  These  artihcial  velocities  clearly  change  the  shape  of  the  large 
vortex  surrounding  the  vane  and  affect  the  radiometric  force  prediction. 


0.20 

0.10 

0.05 

0.00 


Figure  1:  Flow  velocity  and  streamlines  around  a  radiometer  vane  in  a  6  Pa  argon  gas, 
obtained  by  a  finite  volume  approach  to  the  BGK  equation. 


The  third  potential  problem  is  the  lack  of  mass,  momentum  and  energy 
conservation  in  implicit  hnite  volume  schemes  and  the  need  to  use  a  special 
technique  to  impose  conservation  in  explicit  schemes.  The  issues  related  to 
the  conservation  laws  may  be  amplihed  in  transient  hows  for  which  accurate 
prediction  of  time-dependent  properties  is  a  necessity. 

One  possible  way  to  address  the  above  problems  is  to  use  a  DG  type  of 
discretization  both  in  physical  and  velocity  space.  High  order  DG  schemes 
are  generally  formulated  both  for  regular  and  irregular  partitions.  They  have 
compact  stencils  and  are  well  suited  for  adaptive  grid  rehnement.  A  weak 
enforcement  of  huxes  in  DG  discretizations  is  well  suited  for  implementing 
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realistic  gas-surface  interaction  models.  In  particular,  it  does  not  require  the 
use  of  virtual  cells  behind  solid  interfaces.  This  type  of  discretization  will 
generally  allow  one  to  drastically  increase  the  order  of  discretization  and  thus 
reduce  computational  time  in  multi-dimensional  problems.  A  DG  approach 
is  therefore  used  in  this  work. 

Recently  an  application  of  DG  discretizations  to  the  Boltzmann  equation 
was  discussed  in  [21,  18].  In  [19],  the  ellipsoidal-statistical  BGK  equation 
is  solved  using  a  DG  discretization  in  the  spatial  variable  and  discrete  ordi¬ 
nates  in  the  velocity  variable.  In  this  paper  high  order  DG  discretizations  in 
the  spatial  variable  are  augmented  by  high  order  DG  discretizations  in  the 
velocity  variable. 

A  high  order  DG  discretization  in  the  velocity  variable  is  expected  to 
provide  sufficient  accuracy  to  enforce  conservation  of  mass,  momentum  and 
energy  with  no  additional  effort.  It  is  also  important  to  compare  solutions 
obtained  by  high  order  DG  velocity  discretizations  with  the  conservative 
discretization  [12].  Both  steady-state  and  transient  flow  solutions  are  of 
interest. 

3.  The  discontinuous  Galerkin  discretization 

Discontinuous  Galerkin  methods  have  been  applied  to  different  gas  dy¬ 
namic  equations  in  the  past.  A  theoretical  study  of  the  convergence  and 
stability  of  DG  formulations  for  the  Boltzmann  equation  was  conducted  in 
[17].  A  general  description  and  theory  of  Runge-Kutta  DG  methods  may  be 
found  in  [23,  22].  The  first  application  of  a  DG  discretization  of  the  BGK 
equation  in  both  spatial  and  velocity  variables  was  presented  in  [20].  As  in 
[20]  only  approximations  with  hnite  support  are  used  in  this  work.  Let  us 
consider  a  finite  region  T  C  in  the  velocity  space  and  assume  that  the 
distribution  function  f{t,x,u)  and  its  integrals  are  small  in  \  T.  Note 
that  the  method  can  be  extended  to  inhnite  domains  by  using  a  mixture  of 
hnite  and  inhnite  elements  with  appropriate  basis  functions. 

Let  us  introduce  partitions  of  the  spatial  domain  D  G  into  polyhedral 
elements  Kj,  j  =  1, ...  ,N  and  velocity  domain  T  into  polyhedral  elements 
Vi,  i  =  1, . . . ,  M .  Let  polynomial  basis  functions  ippj{x),  p  =  1, . . .  ,k  and 
\i,i{u),  I  =  1, ...  ,s  be  dehned  on  Kj  and  Vi,  respectively.  On  each  phase 
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element  Kj  x  Vi  the  solution  is  sought  in  the  form 

s  At(0 
l=l  p=l 

If  \i^i{u)  are  polynomials  of  degree  not  greater  than  s  and  (Ppj{x)  are  polyno¬ 
mials  of  degree  not  greater  than  k,  then  an  efficient  discretization  is  obtained 
by  setting  /i(/)  =  min(max_degree  —  l,k),  where  max.degree  is  the  desired 
maximum  degree  of  (4).  In  this  paper  max_degree  =  max(A;,  s). 

Let  us  substitute  (4)  into  (1)  and  multiply  the  result  by  a  basis  function 
m  =  1, . . .  ,s.  The  DG  discrete  velocity  formulation  is  obtained  by 
successive  integration  over  Vi  as  follows: 


dtfij{t,  x)  +  d^f'llij{t,  x)  +  dyT1iij{t,  x)  +  % (t,  x) 


=  z/  ID 


'V, 


foXm,i 


(5) 


Here  iij{t,x)  :=  fi,i-pjit)Tp,j{^)  is  the  vector  of  coefficients  of  the  DG 

discrete  velocity  representation  and  square  matrices  T“  = 

(Dj)“^T"  and  =  (Dj)~^TJ’  represent  the  coefficients  of  the  DG  discrete 
velocity  formulation.  Matrix  multiplication  notation  is  used  in  (5)  so  that 

Also, 


TZ,i  =  /  u\Arn,i,  T^l,i  =  /  v\i\m,i, 

JV,  JVi 

'^ml  i  /  DytiI^i  /  Xl^iXjn,i- 

JVi  JVi 


(6) 


In  order  to  derive  the  spatial  discretization,  (6)  is  multiplied  by  a  basis  func- 
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tion  (pgj(x),  q  =  1, . . . ,  /i(m)  and  integrated  over  Kj,  resulting  in 


/i(m}  s  At(0 

p=l  1=1  P=1 

s  m(0  s  pH) 

-  E  U,  E  f'.v.iWC’i,.,  -  E  U,  E  f‘.vjwcki 

1=1  p=i  1=1  p=i 

s  pH)  s  pH) 

+  E  E  + E  E  f-.-ppmCI, 

1=1  p=i  1=1  p=i 


V 


IdKi 


s  „  p(m) 

[(E(^  ml^i  )  ^  -E 

Z  =  1  p=l 


where 


c'  —  r 

^3  ^pq,3 


(Ppj{x)(Pgj{x), 


I  Ki 


f~-iX  _ 

~  ^PQJ 


Tp,j  (2^) ) 


'Ki 


=  ^m,3  =  / 

JKi 


r^z  _ 

^3  ^pg,3 


(Pp,jix)d,^pgj(x), 


'K, 


^dK  _  ^dK 

^3  ~  '-^pgJ 


IdKi 


<fp,j{x)ipgj{x)d(T, 


f^dK  _  z-idK 

^3*  ~  '~"pqJ* 


IdKi 


ippj*{x)ipgj{x)da. 


(8) 


An  upwind  flux  is  used  in  this  formulation.  Thus  j*  is  the  index  correspond¬ 
ing  to  the  element  Kj*  that  shares  the  particular  part  of  the  boundary  with 
Kj.  Projection  operators  Cf^^j  and  appearing  in  (7)  separate  modes 

in  lij  that  propagate  inside  and  outside  Kj,  respectively.  They  are  defined  as 
follows.  Let  A;  and  =  1, . . . ,  s  be  the  eigenvalues  and  the  corresponding 
eigenvectors  of  Tj^-,  respectively.  Let  A^-  be  the  matrix  whose  l-th  column  is 
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We  define  the  diagonal  matrices 

A“  =  diag  (min(0,  Ai), . . . ,  min(0,  A^)), 

A"*"  =  diag  (max(0,  Ai), . . . ,  max(0,  A^)). 

Then  the  projection  operators  are  given  by 

/>-_  A,,A-A-1  />+_  A,,A+A-1 

A  one-dimensional  DG  code  has  been  developed  that  incorporates  the 
above  scheme.  The  application  of  the  code  to  the  shock  wave  and  heat 
transfer  problems  is  discussed  in  the  following  sections. 

4.  The  normal  shock  wave  problem 

The  accuracy  and  convergence  of  the  DG  approach  (8)  is  examined  below 
in  application  to  the  solution  of  the  normal  shock  wave  problem.  This  is 
a  classical  problem  of  gas  dynamics:  as  such  it  was  used  for  accuracy  and 
applicability  analysis  of  a  large  number  of  numerical  methods  of  fluid  dy¬ 
namics.  The  normal  shock  wave  problem  is  characterized  by  a  signihcantly 
non-equilibrium  velocity  distribution  function.  It  allows  one  to  avoid  prob¬ 
lems  associated  with  boundary  conditions  at  solid  walls.  In  this  paper,  a 
weak  shock  wave  in  argon  gas  was  modeled  for  a  Mach  number  of  1.2  and 
an  upstream  number  density  and  temperature  of  1.6x10^^  molecule/m^  and 
300  K,  respectively.  The  argon  viscosity  was  assumed  to  be  2.117“^  N-s/m^ 
at  273  K.  The  length  of  the  computational  domain  was  0.2  m,  which  amounts 
to  about  300  upstream  mean  free  paths.  This  is  large  enough  to  avoid  the 
impact  of  the  size  of  the  computational  domain. 

The  results  obtained  by  the  DG  approach  are  compared  with  the  solution 
of  the  BGK  equation  obtained  by  the  hnite  volume  approach.  (Analysis  of 
the  accuracy  of  the  BGK  approximation  to  the  Boltzmann  equation  for  the 
normal  shock  wave  problem  and  comparison  with  other  kinetic  approaches  is 
not  the  topic  of  the  present  work  and  may  be  found  elsewhere  [15].)  A  hnite 
volume  2D/axisymmetric  code  SMOKE  [16]  has  been  used  to  solve  the  BGK 
equation.  SMOKE  is  a  parallel  code  based  on  conservative  numerical  schemes 
developed  by  L.  Mieussens  [12].  A  second  order  spatial  discretization  is  used 
along  with  explicit  time  integration.  One  thousand  spatial  cells  were  used, 
with  specular  boundaries  set  to  model  ID  how  in  the  longitudinal  direction. 
The  convergence  study  on  the  velocity  grid  was  conducted,  with  the  number 
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of  {x,y,z)  velocity  points  ranging  from  (20,10,10)  to  (60,40,40).  The  results 
for  the  (30,40,40)  grid  are  presented  below. 

Figure  2  shows  the  comparison  of  the  solutions  obtained  by  two  different 
approaches  to  the  model  kinetic  equation.  The  DG  solution  was  obtained 
by  fourth  degree  polynomial  approximations  in  space  on  64  cells  and  eighth 
degree  polynomial  approximations  in  velocity  space  on  32  velocity  bins.  A 
Runge-Kutta  explicit  scheme  was  used  for  the  time  integration.  Only  a  part 
of  the  computational  domain  is  presented  in  order  to  show  more  detail  of  the 
shock  wave  structure.  The  normalized  macroparameters,  U  =  and 

^  +  oo  O  —  CO 

^  rjp _ rjp 

T  =  Tfq - are  shown.  Hereafter  the  DG  solution  is  denoted  by  BGK-DG 

and  the  hnite  volume  solution  is  denoted  by  BGK-FV.  The  results  show  very 
good  agreement  between  the  two  approaches  which  may  be  considered  to  be 
a  verihcation  of  the  developed  DG  code. 


Figure  2:  Normalized  velocity  and  temperature  in  a  Mach  1.2  shock  wave.  Comparison 
between  the  DG  and  FV  approaches. 


High  order  convergence  of  the  DG  method  is  illustrated  in  Figures  3  and  4. 
Figure  3  shows  the  convergence  of  density  prohles  with  respect  to  the  velocity 
variable  for  different  orders  of  the  DG  velocity  approximation.  Hereafter,  the 
relative  error  is  calculated  as  \p  —  pref\/p,  where  p  is  the  gas  density  in  the 
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current  solution  and  Pref  is  the  gas  density  in  a  reference  solution.  Solutions 
in  Figure  3(a)  are  computed  using  DG  approximations  by  first  degree  poly¬ 
nomials,  Figure  3(b)  by  fourth  degree  polynomials  and  those  on  Figure  3(c) 
by  eighth  degree  polynomials.  The  DG  discretization  in  x  is  the  same  in  each 
case  and  used  fourth  degree  polynomials  on  32  cells.  Integration  in  time  was 
conducted  by  a  fifth  order  Runge-Kutta  method.  Note  that  using  different 
numbers  of  spatial  and  velocity  cells  resulted  in  sightly  different  steady-state 
location  of  the  shock  wave  prohle.  The  graphs  were  therefore  obtained  by 
translating  density  prohles  that  correspond  to  different  cell  numbers  so  that 
they  coincided  in  the  center  of  the  shock  wave  where  the  normalized  density 
was  equal  to  0.5.  The  accuracy  of  translation  was  maintained  to  the  eleventh 
digit.  This  explains  the  sharp  drop  in  the  numerical  error  in  the  center  of 
the  shock  wave.  Note  that  the  reference  solutions  pref  in  Figures  3(a)-(c)  are 
different  in  each  case  and  were  obtained  using  the  same  order  of  velocity  dis¬ 
cretization  as  the  compared  solutions  but  on  hner  velocity  grids.  Specifically, 
first  degree  polynomials  on  512  velocity  bins  were  used  as  reference  values  in 
Figure  3(a),  fourth  degree  on  128  bins  in  Figure  3(b)  and  eighth  degree  on 
64  bins  in  Figure  3(c).  This  allows  one  to  show  that  the  solutions  converge 
with  respect  to  the  resolution  in  u  in  each  case.  Moreover,  the  graphs  do 
not  change  signihcantly  if  the  reference  solutions  in  Figures  3(a)  and  (b)  are 
replaced  by  the  solution  obtained  by  the  eighth  degree  polynomials  on  64 
bins  used  in  Figure  3(c).  It  may  therefore  be  concluded  that  the  solutions 
do  converge  to  an  approximate  solution  corresponding  to  a  given  discretiza¬ 
tion  in  X.  At  the  same  time,  comparison  to  a  solution  obtained  on  a  hner 
grid  in  x — for  example,  64  spatial  cells — shows  loss  of  convergence  below  the 
level  of  10“®  at  the  center  of  the  shock  wave,  as  illustrated  in  Figure  3(d). 
This  suggests  that  the  numerical  errors  in  these  solutions  are  dominated  by 
the  errors  of  discretization  in  the  spatial  variable.  In  particular,  no  signih- 
cant  improvement  in  the  solution  can  be  achieved  by  further  velocity  mesh 
rehnement. 

The  shock  wave  solutions  examined  in  Figure  3  exhibit  very  fast  con¬ 
vergence  above  the  10“®  level.  Indeed,  Figures  3(a)  and  (b)  show  an  im¬ 
provement  in  the  solution  accuracy  that  is  seemingly  better  than  the  order 
of  the  polynomial  approximation  used.  This  convergence  can  be  attributed 
to  the  properties  of  the  distribution  function.  It  is  known  that  the  error  of 
evaluating  the  integral  du  by  the  trapezoidal  rule  on  a  sufficiently 

large  hnite  interval  decreases  faster  than  any  power  of  the  mesh  size.  Simi¬ 
larly  fast  convergence  of  quadrature  rules  is  expected  for  all  smooth  functions 
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that  rapidly  decrease  at  infinity.  It  is  believed  that  the  fast  convergence  of 
the  macroparameters  with  respect  to  the  resolution  in  velocity  space  is  due  to 
the  exponential  decrease  of  the  distribution  function.  The  fast  convergence 
is  not  limited  to  high  order  DG  approximations.  It  is  expected  to  hold  for 
the  discrete  ordinate  approach  on  uniform  velocity  bins  as  well.  The  other 
necessary  condition  for  fast  convergence  is  the  smoothness  of  the  solution. 
In  particular,  one  should  not  expect  fast  convergence  in  velocity  space  if  the 
solution  is  contaminated  by  errors  of  discretization  in  physical  space.  Note 
that  the  convergence  becomes  only  second  order  in  Figure  3(b)  once  the 
numerical  error  had  reached  the  level  of  the  error  of  spatial  discretization. 
Interestingly,  eighth  degree  polynomial  approximations  are  much  less  sensi¬ 
tive  to  the  roughness  of  the  solution  and  converge  to  almost  the  round-off 
level. 

The  CPU  times  that  correspond  to  the  simulations  shown  in  Figure  3 
are  given  in  Table  1.  The  results  were  obtained  on  an  AMD  Opteron  252 
processor.  Analysis  of  the  Table  1  and  Figure  3  shows  that  it  is  impractical 
to  use  more  than  64  bins  in  the  first  degree,  32  bins  in  the  fourth  degree  and 
16  bins  in  the  eighth  degree  polynomial  approximations.  Thus  the  required 
CPU  times  are  13,060  seconds  for  s  =  1,  15,486  seconds  for  s  =  4,  and 
15,  757  seconds  for  s  =  8.  While  these  run  times  are  comparable,  the  s  =  8 
solution  may  be  preferred  since  it  has  a  faster  convergence. 


s  =  1 

s  =  4 

s  =  8 

M 

t,  CPU  sec. 

M 

t,  CPU  sec. 

M 

t,  CPU  sec. 

32 

6,349 

8 

3,891 

8 

7,713 

64 

13,060 

16 

7,683 

16 

15,757 

128 

25,232 

32 

15,486 

32 

30,117 

256 

49,665 

64 

30,579 

64 

61,188 

Table  1:  CPU  time  as  a  function  of  degree  of  polynomial  approximation  in  the  velocity 
space. 


Figure  4  illustrates  the  convergence  of  the  density  profile  with  respect 
to  the  spatial  variable.  Two  cases  of  DG  approximation  in  physical  space 
are  shown,  with  first  degree  polynomials  plotted  in  Figure  4(a)  and  fourth 
degree  polynomials  in  Figure  4(b).  The  temporal  Runge-Kutta  integration 
is  second  and  fifth  order,  respectively.  The  DG  velocity  discretization  is 
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Figure  3:  Convergence  of  the  DG  solution  in  a  Mach  1.2  shock  wave  as  a  function  of  the 
number  of  velocity  bins.  Comparison  of  the  DG  discretizations  by  (a)  first,  (b)  fourth  and 
(c)  eighth  degree  polynomials.  In  figure  (d)  solutions  on  32  spatial  cells  are  compared  to 
a  solution  on  64  spatial  cells. 


conducted  using  eighth  degree  polynomials  on  32  uniform  velocity  bins  in 
each  case.  Note  that  the  scheme  using  fourth  order  polynomials  exhibits 
considerably  less  oscillation  and  has  better  convergence  at  the  center  of  the 
shock  wave.  In  both  cases,  the  convergence  generally  corresponds  to  the 
order  of  the  polynomial  approximation  used. 

The  presented  results  show  that  the  DG  discretization  produces  accurate 
solutions  that  converge  with  high  order  both  in  physical  and  velocity  vari¬ 
ables.  Note  also  that  since  high  order  discretization  in  velocity  space  does 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


4  THE  NORMAL  SHOCK  WAVE  PROBLEM 


Figure  4:  Convergence  of  density  in  a  Mach  1.2  shock  wave  as  a  function  of  spatial  cells. 
Comparison  of  (a)  second  and  (b)  fifth  order  DC  discretizations. 


not  lower  the  CFL  condition  for  the  temporal  discretization,  it  may  serve  as 
an  inexpensive  way  to  increase  solution  accuracy. 

Implementation  of  the  DG  procedure  for  the  velocity  and  physical  space 
discretizations  is  more  cumbersome  than  that  of  a  typical  FV  approach.  As 
a  result,  computational  time  for  the  DG  method  is  larger  than  that  of  FV 
methods  with  the  same  order  of  discretizations.  In  the  case  of  the  second 
order  in  space  and  hrst  order  in  velocity  discretizations  the  FV  solver  SMOKE 
was  two  to  three  times  faster  than  the  DG  code.  Performance  of  the  DG 
method,  however,  can  be  improved  by  increasing  the  order  of  polynomial 
approximation. 

Table  2  shows  the  maximum  error  near  the  center  of  the  shock  wave  and 
the  GPU  time  used  for  simulations  shown  in  Figure  4.  It  can  be  seen  that  the 
simulations  that  use  a  piece-wise  linear  approximation  in  space  converge  to 
~  10“^  for  128  cells.  At  the  same  time,  similar  error  is  achieved  for  only  16 
cells  when  a  fourth  degree  polynomial  approximation  is  used.  Gomparison  of 
computational  times  shows  that  the  high  order  method  is  about  3.77  times 
faster.  Table  2  also  suggests  that  more  accurate  simulations  will  require  con¬ 
siderably  longer  time  for  low  order  techniques.  At  least  four  additional  grid 
rehnements  are  necessary  for  the  piece- wise  linear  simulations  to  achieve  10“® 
accuracy.  Since  each  rehnement  results  in  about  3.5  increase  in  GPU  time. 
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simulations  based  on  first  degree  polynomial  approximations  will  be  at  least 
150  times  slower  than  those  for  fourth  degree  approximations.  Performance 


k  = 

1 

k  = 

4 

N 

rel.  err. 

t,  GPU  sec. 

N 

rel.  err. 

t,  GPU  sec. 

16 

8.2E-3 

985.97 

4 

1.3E-2 

876.09 

32 

4.1E-3 

3264.48 

8 

3.8E-3 

3151.16 

64 

l.OE-3 

11551.34 

16 

4.1E-4 

11825.48 

128 

1.2E-4 

44593.05 

32 

2.2E-6 

45416.36 

Table  2:  CPU  times  for  simulations  by  second  order  and  fifth  order  methods  shown  on 
Figure  4. 


of  the  DG  method  (7)  may  be  further  improved  by  careful  selection  of  ba¬ 
sis  functions  and  quadrature  formulas  for  the  evaluation  of  moments  of  the 
distribution  function.  One  way  to  do  that  is  to  generalize  the  approach  of 
[21]  to  the  DG  methods  and  use  Lagrange  polynomials  on  Hermite  nodes  as 
velocity  basis  functions.  The  advantage  of  such  a  basis  is  that  matrices  (6) 
are  diagonal.  This  will  simplify  implementation  and  reduce  computational 
time. 

5.  Heat  transfer  between  parallel  plates 

The  normal  shock  wave  problem,  while  being  important  for  testing  the 
method  applicability  and  convergence  for  gas  flow  modeling,  does  not  allow 
one  to  analyze  the  approach  accuracy  and  robustness  for  modeling  gas-solid 
interfaces.  However,  gas-surface  interactions  are  extremely  important  for 
low-speed  gas  flows  in  microscale  devices  where  reliable  modeling  of  such  in¬ 
teractions  is  indispensable.  The  Maxwell  model  of  gas-surface  collisions  and, 
in  particular,  collisions  with  full  momentum  and  energy  accommodation  at 
the  surface  (fully  diffuse  reflection)  are  most  widely  used  in  computations 
of  microscale  gas  flows  at  the  kinetic  level.  The  Maxwell  model  is  naturally 
suitable  for  particle  approaches  such  as  DSMG.  Deterministic  kinetic  meth¬ 
ods,  both  hnite  difference  and  hnite  volume,  generally  suffer  from  the  loss  of 
accuracy  related  to  the  discontinuity  in  the  molecular  velocity  distribution 
function  at  the  surface.  An  exception  here  is  the  fully  specular  reflection, 
which  is  not  physically  realistic. 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


5  HEAT  TRANSFER  BETWEEN  PARALLEL  PLATES 


In  order  to  study  the  affect  of  gas-surface  interactions  on  accuracy  and 
convergence  of  the  DG  approach,  a  one- dimensional  heat  transfer  between 
parallel  plates  was  considered,  with  fully  diffuse  reflection  assumed  at  the 
wall.  The  cold  and  hot  wall  temperatures  were  300  K  and  1000  K,  respec¬ 
tively.  Argon  gas  was  modeled  for  two  Knudsen  numbers,  Kn  =  0.01  and 
1. 

First,  we  examine  the  conservation  of  the  total  mass  in  the  numerical 
solution.  In  the  previous  section  it  was  concluded  that  low  order  velocity 
discretization  methods  may  provide  accurate  estimates  of  moments  if  the 
numerical  solution  is  sufficiently  smooth.  While  this  conclusion  is  true  for 
the  shock  wave  problem,  it  is  not  the  case  for  the  heat  transfer  problem. 
Figure  5  shows  the  relative  error  in  the  total  mass  of  the  gas  between  the 
plates  as  a  function  of  time.  Solutions  obtained  by  fourth  and  eighth  degree 
polynomial  approximations  in  velocity  are  shown  in  Figures  5(a)  and  (b), 
respectively.  The  spatial  DG  discretization  is  by  fourth  degree  polynomials 
with  32  spatial  cells  in  both  cases. 

As  expected,  the  error  in  the  total  mass  decreases  as  the  number  of  ve¬ 
locity  bins  increases.  However  the  rate  of  decrease  suggests  that  for  a  hfth 
order  polynomial  approximation  the  decrease  in  error  is  only  of  second  or¬ 
der.  At  the  same  time,  the  eighth  degree  polynomial  approximation  shown 
on  Figure  4(b)  conserves  mass  within  the  ninth  digit  even  when  only  16  ve¬ 
locity  bins  are  used.  Note  that  the  rate  of  convergence  is  lost  below  the 
10“®  level.  However,  at  this  point,  the  simulations  had  reached  the  rounding 
off  error  limit  and  any  further  improvement  is  not  expected.  Note  that  the 
use  of  eighth  degree  polynomial  approximations  also  allows  one  to  conserve 
energy  flux  and  maintain  zero  mass  and  momentum  fluxes  at  the  10“®  level 
of  accuracy  throughout  the  computational  domain. 

Similar  observations  can  be  made  regarding  the  convergence  of  the  density 
prohle  with  respect  to  the  velocity  variable.  Gomparison  of  solutions  obtained 
by  fourth  and  eighth  degree  polynomial  approximations  in  velocity  space  is 
given  in  Figure  6.  The  cold  wall  is  located  at  a;  =  0  and  the  hot  wall  is 
located  at  a;  =  0.1.  As  seen  in  Figure  6(a)  the  solutions  obtained  by  the 
fourth  degree  approximation  converge  with  the  third  order  or  more  slowly. 
As  in  the  mass  conservation  analysis,  the  order  of  convergence  is  two  units 
less  than  the  order  of  polynomial  approximation  used.  It  is  reasonable  to 
assume  that  the  convergence  is  slower  because  the  solution  is  not  sufficiently 
smooth  (the  reasons  for  this  are  discussed  below).  At  the  same  time,  the 
results  obtained  using  eighth  degree  polynomials  shown  in  Figure  6(b)  are 
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Figure  5:  Convergence  of  relative  error  in  the  total  mass  with  respect  to  the  resolution  in 
velocity  space.  Comparison  of  (a)  fourth  and  (b)  eighth  degree  polynomial  approximations. 


not  sensitive  to  the  non-smoothness  of  the  solution  and  quickly  converge  to 
the  level  of  the  round-off  error. 


Figure  6:  Convergence  of  the  solution  in  a  Kn=0.01  heat  transfer  problem  as  a  function  of 
the  number  of  velocity  bins.  Comparison  of  (a)  fourth  and  (b)  eighth  degree  polynomial 
DC  approximation  in  velocity  space. 


It  is  not  immediately  obvious  why  the  solution  looses  its  smoothness  in 
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the  entire  domain.  A  study  of  convergence  with  respect  to  the  resolution 
in  physical  space  suggests  that  the  non-smoothness  is  caused  by  the  diffuse 
boundary  conditions.  The  convergence  of  the  density  prohle  with  respect  to 
the  resolution  in  physical  space  is  shown  in  Figure  7,  where  solutions  obtained 
using  hrst  and  fourth  degree  polynomial  DG  approximations  in  x  are  shown. 
The  Runge-Kutta  time  integrations  are  second  and  hfth  order,  respectively. 
The  discretization  in  velocity  space  is  conducted  on  32  velocity  bins  using 
eighth  degree  polynomials  for  both  cases.  The  order  of  convergence  for  the 
density  prohle  is  second  order  at  best.  Also,  there  is  an  oscillation  near  the 
walls.  We  believe  that  the  reason  for  this  is  the  discontinuity  in  the  velocity 
distribution  function  near  the  wall,  as  is  explained  below. 

Generally  the  solution  to  the  heat  transfer  problem  at  a  hnite  Knudsen 
number  exhibits  a  temperature  jump  at  the  wall.  This  means  that  the  tem¬ 
perature  of  gas  molecules  that  collide  with  the  wall  and  the  temperature  of 
the  wall  are  not  the  same.  Noticing  that  in  the  diffuse  boundary  conditions 
the  rehected  gas  molecules  are  modeled  by  the  Maxwellian  distribution  func¬ 
tion  at  the  wall  temperature  and  the  colliding  molecules  are  distributed  very 
close  to  a  Maxwellian  but  have  a  different  temperature,  the  diffusion  bound¬ 
ary  conditions  force  the  solution  to  form  a  discontinuity  at  the  wall  at  m  =  0. 
Thus  high  gradients  will  appear  in  the  direction  of  the  velocity  variable  near 
the  wall  at  point  u  =  0.  These  gradients  propagate  inside  the  domain  and 
decay  quickly  with  a  rate  governed  by  the  collision  frequency  u{x,t).  This 
in  turns  creates  high  gradients  in  the  x  direction,  which  produce  a  high  fre¬ 
quency  error  in  the  solution.  The  fact  that  the  cold  boundary  has  a  higher 
density  and  thus  effectively  a  higher  collision  frequency  explains  why  more 
oscillations  are  observed  near  the  cold  wall  than  near  the  hot  wall. 

Let  us  now  compare  the  macroparameters  obtained  by  the  FV  and  DG 
methods  of  discretization.  The  ratios  of  the  FV  temperature  prohles  to  the 
corresponding  reference  DG  solution  computed  using  the  ninth  order  in  ve¬ 
locity  and  hfth  order  in  physical  space  scheme  are  presented  in  Figure  8(a)  for 
Kn=l.  The  hgure  shows  that  the  FV  solution  has  converged  with  respect  to 
the  number  of  spatial  cells,  since  the  solution  does  not  change  when  the  num¬ 
ber  of  cells  is  increased  from  100  to  300.  However,  there  still  is  a  difference 
between  the  FV  and  DG  results  that  is  related  to  the  hnite  number  of  veloc¬ 
ity  bins  used  in  the  FV  scheme  (in  this  case,  (40,30,30)).  The  FV  scheme 
corresponds  to  the  discrete  ordinate  approach  on  a  uniform  velocity  grid  and, 
as  in  low  order  DG  schemes  (see  Figure  6),  requires  a  very  large  number  of 
velocity  bins  to  achieve  convergence  better  than  in  the  second  digit.  The 
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(b) 


Figure  7:  Convergence  of  density  in  a  Kn=0.01  heat  transfer  problem  as  a  function  of 
the  number  of  cells.  Comparison  of  (a)  first  degree  and  (b)  fourth  degree  polynomial  DC 
discretizations. 


main  reason  for  the  inability  of  the  FV  method  to  accurately  capture  the 
temperature  prohle  for  a  limited  number  of  velocity  bins  is  primarily  related 
to  the  discontinuity  of  the  velocity  distribution  function  similar  to  that  men¬ 
tioned  earlier.  For  the  relatively  high  Knudsen  number  that  we  have  used, 
the  distribution  function  at  any  given  spatial  location  represents  a  combi¬ 
nation  of  two  half-Maxwellians  that  correspond  to  particles  reflected  on  the 
cold  and  hot  plates.  A  large  number  of  velocity  bins  is  needed  in  the  FV 
method  in  order  to  properly  capture  such  a  discontinuous  shape. 

Another  important  metric  of  accuracy  of  the  obtained  solution  is  the 
flow  velocity  between  the  plates.  In  the  absence  of  plate  motion,  the  gas 
bulk  velocity  should  be  zero  at  any  point  between  the  plates.  The  values  of 
average  gas  velocity  in  the  direction  perpendicular  to  the  plates  are  given 
in  Figure  8(b)  for  the  FV  and  DG  approaches.  The  FV  approach  that  uses 
100  spatial  cells  is  characterized  by  relatively  high,  on  the  order  0.1  m/s, 
velocities  in  the  regions  near  the  walls.  The  flow  in  the  central  region  also 
has  noticeable  velocities,  with  gas  moving  at  speeds  about  0.002  m/s  from 
the  hot  to  the  cold  wall.  Magnitudes  of  the  flow  velocities  are  relatively 
small  compared  to  the  thermal  velocity  of  about  350  m/s.  However,  these 
magnitudes  may  be  considerably  high  for  some  low-speed  microscale  flow 
applications.  The  increase  in  the  number  of  cells  from  100  to  300  allowed  a 
signihcant  reduction  of  flow  velocities,  although  they  are  still  visible  both  in 
the  central  region  and  especially  near  the  walls.  It  is  interesting  to  note  that 
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this  numerical  effect,  associated  with  the  boundary  condition  specification  in 
the  FV  method,  is  observed  both  at  the  cold  and  at  the  hot  wall;  the  error 
near  the  cold  wall  is  somewhat  higher.  For  the  DG  discretization,  the  error 
in  the  flow  velocity  calculation  is  on  the  order  of  10“®  m/s,  and  therefore 
is  not  visible  in  the  figure  (the  DG  result  appears  to  be  a  straight  line  at 
U,  =  0). 


(b) 


Figure  8:  FV  to  DG  temperature  ratio  (a)  and  gas  flow  velocity  (b),  for  heat  transfer 
problem  at  Kn=l. 


The  decrease  in  the  Knudsen  number  from  1  to  0.01  results  in  a  signifi¬ 
cantly  faster  collisional  relaxation  of  the  distribution  function,  which  is  close 
to  the  Gaussian  shape.  A  smaller  number  of  velocity  bins  is  therefore  needed 
in  the  FV  method  to  accurately  describe  the  flow.  This  is  clearly  seen  in  Fig¬ 
ure  9(a),  where  the  ratios  of  the  FV  to  DG  temperature  profiles  are  shown. 
At  the  same  time,  more  spatial  cells  are  required  to  properly  capture  the  flow 
gradients.  When  100  cells  were  used,  the  error  was  over  one  percent.  This  is 
not  surprising  since  the  cell  size  on  the  order  of  the  gas  mean  free  path  is  not 
expected  to  produce  accurate  results  in  this  low  order  scheme.  For  the  case 
of  1000  cells  the  error  was  significantly  smaller  and  was  less  than  0.5%  even 
near  the  cold  surface  where  the  error  was  largest.  For  further  error  reduction, 
more  velocity  bins  are  needed  than  the  used  value  of  (40,30,30).  However, 
as  was  mentioned  above,  the  accuracy  of  the  FV  approach  is  limited  by  the 
approximations  used  for  the  boundary  conditions  applied  at  the  wall. 

The  issue  with  the  boundary  conditions  is  illustrated  in  Figure  9(b).  For 
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the  DG  discretization,  the  flow  velocity  is  fairly  close  to  zero,  although  the 
error  is  somewhat  larger  than  for  Kn  =1.  As  discussed  earlier,  the  most 
noticeable  error  for  this  method  is  observed  near  the  cold  surface,  where  the 
magnitude  of  the  flow  velocity  amounts  to  about  0.005  m/s.  In  the  rest  of 
the  computational  domain  it  is  on  the  order  of  10“^  m/s.  The  error  was 
considerably  larger  for  the  FV  method,  even  when  1000  cells  were  used;  in 
this  case  the  cell  size  was  about  10%  of  the  gas  mean  free  path.  Note  that  the 
error  in  flow  velocity  was  not  localized  near  the  walls,  but  propagated  into 
the  central  region  of  the  domain,  where  it  was  over  two  orders  of  magnitude 
larger  than  in  the  DG  approach.  The  error  was  even  larger  for  100  cells, 
where  it  reached  values  in  excess  of  1  m/s. 


(b) 


Figure  9:  FV  to  DG  temperature  ratio  (a)  and  gas  flow  velocity  (b),  for  heat  transfer 
problem  at  Kn=0.01. 


Analysis  of  transient  flow  evolution  for  the  heat  transfer  problem  was  also 
conducted.  At  several  time  moments,  the  FV  and  DG  results  were  compared 
for  the  two  Knudsen  numbers  under  considerations.  The  results  were  found 
to  agree,  with  an  error  close  to  that  observed  for  the  steady  state  calculations 
presented  earlier. 

6.  Conclusions 

High  order  discontinuous  Galerkin  discretizations  both  in  spatial  and  ve¬ 
locity  variables  were  applied  to  the  BGK  equation.  The  normal  shock  wave 
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and  the  heat  transfer  between  parallel  plates  were  analyzed  in  detail.  It  was 
found  that  the  solutions  to  the  normal  shock  wave  problem  exhibit  very  fast 
convergence  with  respect  to  the  resolution  in  the  velocity  variable.  The  con¬ 
vergence  of  the  solutions  with  respect  to  spatial  variable  corresponds  to  the 
order  of  polynomial  interpolation.  Furthermore,  it  was  found  that  high  order 
DG  discretizations  of  the  kinetic  equations  conserve  mass,  momentum  and 
energy  to  a  high  precision.  The  DG  solutions  were  compared  to  solutions 
obtained  by  a  FV  conservative  technique  and  were  found  to  be  in  excellent 
agreement.  High  order  discretizations  in  the  velocity  variable,  such  as  those 
using  basis  polynomials  of  eighth  degree,  were  found  to  produce  results  not 
sensitive  to  the  roughness  in  the  solution. 

It  has  been  observed,  however,  that  the  solutions  to  the  heat  transfer 
problem  do  not  exhibit  a  proper  convergence  rate  with  respect  to  the  spatial 
variable.  The  solutions  computed  using  approximations  by  polynomials  of 
first  and  fourth  degrees  converge  with  the  second  order  or  more  slowly.  Loss  of 
convergence  in  the  heat  transfer  problem  is  attributed  to  the  diffuse  boundary 
conditions.  It  appears  that  a  standard  formulation  of  the  DG  method  is  not 
well  suited  to  handle  the  discontinuity  at  the  boundary  that  is  intrinsic  to 
diffuse  boundary  conditions.  Development  of  an  improved  DG  techniques 
that  handles  such  a  discontinuity  is  therefore  necessary. 
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